Confidence in Education

---

This document aims to explores the relationship between confidence in the US education system and actual government spending on education.

Loading Packages and Data

In this section, we pull in data downloaded from the Bureau of Economic Analysis (BEA) website and from the General Social Survey (GSS).

### Load packages
library(tidyverse)
library(readxl)
require(forecast)
library(zoo)
library(kableExtra)
library(astsa)

### Load data
gved_y_raw <- read_excel("data/nipa_3155.xlsx")
gv_y_raw   <- read_excel("data/nipa_116a.xlsx")
gv_q_raw   <- read_excel("data/nipa_116q.xlsx")
pop_raw    <- read_excel("data/nipa_21q.xlsx")
gss_raw    <- read_excel("data/gss_educ.xlsx")

Cleaning data

First, lets clean the data on government spending over time in the US. Unfortunately, the BEA’s data on inflation adjusted spending on education only goes back to 2007. To work around this, we can take the ratio of nominal education spending to nominal total government spending and multiply real total government spending to get an estimate of real total government spending on education. However, total government spending on education is still misleading since it is increasing with things like population growth. To fix this, we can calculate education spending per capita by dividing by the U.S. population. We create two dataframes for the analysis: one with yearly data that we can use to compare to the yearly GSS data, and another that is quarterly for time series analysis.

### BEA NIPA data 
## Population to get per capita spending 
pop <- pop_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,2,45)) %>%
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, quarter = V2, popq = V3) %>%
  fill(year, .direction = "down") %>% 
  mutate(yearq = paste0(year,quarter),
         year  = as.numeric(year),
         popq = as.numeric(popq)) %>%
  select(year, yearq, popq)

pop_y <- pop %>%
  group_by(year) %>%
  summarise(popy = mean(popq, na.rm = TRUE))

## Getting ratio of education spending to total
govt_educ_y <- gved_y_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,3,31)) %>%
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, govt = V2, educ = V3) %>%
  mutate(educ = as.numeric(educ),
         govt = as.numeric(govt),
         year = as.numeric(year)) %>%
  mutate(educ_govt = educ / govt)

## Getting annual government spending 
govt_y <- gv_y_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,24)) %>% 
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, govty = V2) %>%
  mutate(year = as.numeric(year),
         govty = as.numeric(govty)) %>%
  left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
  filter(!is.na(educ_govt)) %>%
  mutate(gved_y = as.numeric(govty) * as.numeric(educ_govt)) %>%
  inner_join(pop_y, by = "year") %>%
  mutate(gved_pc_y = (gved_y * 1e6) / (popy * 1e3)) %>%
  select(year, gved_pc_y)
  
## Quarterly government spending on education (for time series forecast)
govt_q <- gv_q_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,2,24)) %>% 
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, quarter = V2, govtq = V3) %>%
  fill(year, .direction = "down") %>% 
  mutate(yearq   = paste0(year,quarter),
         year    = as.numeric(year),
         quarter = as.numeric(substr(quarter,2,2)),
         govtq   = as.numeric(govtq)) %>%
  select(c(year,quarter,yearq,govtq)) %>%
  left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
  filter(!is.na(educ_govt)) %>%
  mutate(gved_q = govtq * educ_govt) %>%
  inner_join(pop, by = c("yearq","year")) %>%
  mutate(gved_pc_q = (gved_q * 1e6) / (popq * 1e3),
         yearq = as.yearqtr(yearq)) %>%
  select(c(year,quarter,yearq, gved_pc_q))

Regarding the variables from the GSS, we use the following surveys:

• The coneduc survey asks whether the respondent has a great deal of confidence, only some confidence, or hardly any confidence in the institution of education in the US.
• The nateduc survey asks if the respondent thinks the US government is spending too much, too little, or the right amount on education.

We create two dataframes, one that has the average response per year and another that has the number of responses per response per year.

### GSS Education variables
## Taking an average of responses per year
gss_num <- gss_raw %>%
  mutate(nateduc_num = case_when(
            nateduc == "TOO LITTLE" ~ 1,
            nateduc == "ABOUT RIGHT" ~ 2,
            nateduc == "TOO MUCH" ~ 3,
            TRUE ~ NA_real_
          ),
         coneduc_num = case_when(
            coneduc == "HARDLY ANY" ~ 1,
            coneduc == "ONLY SOME" ~ 2, 
            coneduc == "A GREAT DEAL" ~ 3,
            TRUE ~ NA_real_
            ),
         year = as.numeric(year)) %>%
  group_by(year) %>% 
  summarise(mean_nateduc = mean(nateduc_num, na.rm = TRUE),
            mean_coneduc = mean(coneduc_num, na.rm = TRUE),
            n_obs = n()) %>%
  ungroup() %>%
  mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE))

## Counting how many responses of each type per year
gss_count <- gss_raw %>%
  mutate(nateduc = ifelse(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH"), nateduc, NA),
         coneduc = ifelse(coneduc %in% c("A GREAT DEAL","HARDLY ANY","ONLY SOME"), coneduc, NA)) %>%
  select(-id_) %>%
  mutate(year = as.numeric(year))  %>%
  pivot_longer(
    cols = c(nateduc, coneduc),
    names_to = "question",
    values_to = "response"
  ) %>%
  group_by(year, question, response) %>%
  summarise(n = n(), .groups = "drop")

Analysis

First, lets look at how government spending on education per capita has evolved over time.

## NIPA Data
ggplot(data = govt_y, aes(x=year,y=gved_pc_y)) +
  geom_line() + 
  labs(title="Real Government Spending on Education Per Capita",
       x = "Year",
       y = "Chained (2017) Dollars") +
  theme_minimal()

This chart shows that the government has been spending more per capita on education over time. There were some dips in the 1980s and in the early 2010s, but overall it has been increasing. However, spending has stagnated a bit since the early 2000s.

Now that we’ve seen how spending has changed over time, lets look at confidence in education. We can make line charts for the proportion of each response the coneduc question.

## GSS Data
gss_prop <- gss_count %>%
  filter(!is.na(response)) %>%
  group_by(year, question) %>%
  mutate(prop = n / sum(n)) 

# Confidence in Education
coneduc_prop <- gss_prop %>%
  filter(question == "coneduc")

ggplot(coneduc_prop, aes(x = year, y = prop, color = response)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "How much confidence do you have in the education\nsystem in the US?",
    x = "Year",
    y = "Proportion of Respondents",
    color = "Response"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5))

It appears that confidence in the education system has decreased over time. The proportion of respondents who have a great deal of confidence has decreased from around 40% to less than 20% and the proportion who have hardly any confidence has increased from less than 10% to around 25%.

We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypotheses are:

\(H_0:\) The distribution of responses to how much confidence in education is the same prior to 1990 vs after 1990.

\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.

# Choose two periods for comparison (example: early vs recent)
gss_test1 <- gss_raw %>%
  filter(coneduc %in% c("HARDLY ANY","ONLY SOME","A GREAT DEAL")) %>%
  mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
  group_by(period, coneduc) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = coneduc, values_from = n)

gss_test1

# A tibble: 2 × 4
  period `A GREAT DEAL` `HARDLY ANY` `ONLY SOME`
  <chr>           <int>        <int>       <int>
1 Early            7076         2366       11313
2 Late             7206         5488       17115

# Run the chi-square test
chisq.test(gss_test1[ , -1])   

    Pearson's Chi-squared test

data:  gss_test1[, -1]
X-squared = 831.81, df = 2, p-value < 2.2e-16

Since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.

What about feelings about education spending? Let’s make a similar line chart over time to explore this.

# Spending on Education
nateduc_prop <- gss_prop %>%
  filter(question == "nateduc")

ggplot(nateduc_prop, aes(x = year, y = prop, color = response)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Is US spending on education too much, too little, or about right?",
    x = "Year",
    y = "Proportion of Respondents",
    color = "Response"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5))

Since the beginning of the survey, most respondents believed that the government was spending too little on education, and this proportion has increased significantly to almost 80% of respondents.

We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypothese are:

\(H_0:\) The distribution of responses to how much are we spending on education is the same prior to 1990 vs after 1990.

\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.

# Choose two periods for comparison (example: early vs recent)
gss_test2 <- gss_raw %>%
  filter(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH")) %>%
  mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
  group_by(period, nateduc) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = nateduc, values_from = n)

gss_test2

# A tibble: 2 × 4
  period `ABOUT RIGHT` `TOO LITTLE` `TOO MUCH`
  <chr>          <int>        <int>      <int>
1 Early           6183        10351       1440
2 Late            4989        17102       1432

# Run the chi-square test
chisq.test(gss_test2[ , -1])   

    Pearson's Chi-squared test

data:  gss_test2[, -1]
X-squared = 1064.8, df = 2, p-value < 2.2e-16

Again, since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.

Finally, lets plot the mean of these two variables over time to create an even more clear visualization of the decline in sentiment.

gss_long <- gss_num %>%
  pivot_longer(
    cols = c(mean_nateduc, mean_coneduc),
    names_to = "variable",
    values_to = "mean_value"
  )

ggplot(gss_long, aes(x = year, y = mean_value, color = variable)) +
  geom_line(size = 1.1) +
  labs(
    title = "Mean Education Attitudes Over Time",
    x = "Year",
    y = "Mean Value",
    color = "Variable"
  ) +
  theme_minimal()

Now that we’ve looked at the variables of interest separately, lets investigate how spending and confidence in education relate to each other. The following charts show how opinions about the amount of education spending relate to actual education spending.

nateduc_nipa <- gss_prop %>%
  filter(question == "nateduc") %>%
  left_join(govt_y, by = "year")

ggplot(nateduc_nipa, aes(x = gved_pc_y, y = prop)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~ response) +
  labs(
    title = "Confidence in amount of Education Spending vs Actual Spending",
    x = "Real Education Spending Per Capita",
    y = "Proportion of Respondents"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

The negative slope for the respondents who believe the government spends about the right amount on education and the positive slope for the those who believe the government spends too little on education indicate that even when the government prioritize education more, people still believe that the government isn’t spending enough. This may indicate that despite some efforts to spend more, the government still isn’t spending as much on education as they should.

Doing the same with confidence in the education system gives similar results.

coneduc_nipa <- gss_prop %>%
  filter(question == "coneduc") %>%
  left_join(govt_y, by = "year")

ggplot(coneduc_nipa, aes(x = gved_pc_y, y = prop)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~ response) +
  labs(
    title = "Confidence in the US Education System vs Spending",
    x = "Real Education Spending Per Capita",
    y = "Proportion of Respondents"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

Finally, we can do the same thing but using the mean of the surveys.

gss_num_govt <- gss_num %>%
  inner_join(govt_y, by = "year") %>%
  slice(-1)

# Mean_nateduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_nateduc)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "National Education Spending Attitudes vs Education Spending",
    x = "Govt Spending per Capita",
    y = "Mean National Education Attitude"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

# Mean_coneduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_coneduc)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Confidence in Education vs Education Spending",
    x = "Govt Spending per Capita",
    y = "Mean Confidence in Education"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

All these plots serve to show that even as government spending on education increases, people still increasingly have lost trust in the education system in the US. Let’s explore how these trends may evolve in the future with time series forecasts.

ARIMA

In this section, we analyze and attempt to forecast the time series for education spending and confidence in education. We focus only on confidence in education rather than feelings about speding on education since it more closely relates to our overall topic of trust in american institutions.

Government Spending on Education

gv_ts <- ts(govt_q$gved_pc_q,
            start = c(min(govt_q$year),min(govt_q$quarter)),
            frequency = 4) 

Lets start with a classical decomposition.

gv_ts_decomp <- decompose(gv_ts)
autoplot(gv_ts_decomp)

The trend component of the chart shows a clear upward trend over time that matches the data. The seasonal component is also showing clear seasonality with spikes that appear to occur at the beginning of each year, however the effect of the seasonality seems small compared to the trend. Finally, the remainder shows fluctuations around zero that aren’t very uniform. It doesn’t seem likely that the remainder is showing cyclical patterns like recessions since it doesn’t seem to follow any particular pattern, but perhaps it reflects other outside political pressures on spending that are difficult to predict.

Let’s plot the lag and autocorrelation function plots to investigate the patterns more closely.

gglagplot(gv_ts, do.lines=FALSE, lags = 9) +
  xlab("Lags") + ylab("Yt") + 
  ggtitle("Lag Plot for education spending") + 
  theme(axis.text.x=element_text(angle=45, hjust=1))

The linear relationships in the lag plots suggest high serial correlation among consecutive observations.

ggAcf(gv_ts, 50) +
  ggtitle("Autocorrelation Function for Real\nGovernment Spending on Education Per Capita") +
  theme_minimal()

This plot shows clear autocorrleation in the data, however the seasonality isn’t as present. Lets confirm stationarity with an augmented Dickey Fuller test.

tseries::adf.test(gv_ts)

    Augmented Dickey-Fuller Test

data:  gv_ts
Dickey-Fuller = -2.4123, Lag order = 6, p-value = 0.4023
alternative hypothesis: stationary

Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our suspicions from the ACF plot.

We can try to correct this with differencing.

diff_1 <- diff(gv_ts)
ggAcf(diff_1,50) +
  ggtitle("ACF of First-Order Differenced Series") +
  theme_minimal()

It appears that one round of differencing may not have been sufficient since there is still a decreasing pattern every 4 lags. Let’s try one more time to be safe.

diff_2 <- diff(gv_ts, differences = 2)
ggAcf(diff_2,50) +
  ggtitle("ACF of Second-Order Differenced Series") +
  theme_minimal()

ggPacf(diff_2,50) +
  ggtitle("PACF of the Differenced Series") +
  theme_minimal()

This looks a little bit better, perhaps we will need to test for both 1 and 2 lags in the Arima. With this, we can proceed with the ARIMA.

First lets find the optimal parameters. Based on the differenced ACF and PACF plot, I believe the correct range for q is 0-8 and for p is 0-4.

# Define parameter ranges
p_range <- 0:4
d_range <- 1:2
q_range <- 0:9

# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)

# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)

# Initialize index for matrix row
i <- 1

# Loop through combinations of ARIMA model parameters
for (d in d_range) {
  for (q in q_range) {
    for (p in p_range) {
  
      # Fit ARIMA model with specified (p,d,q)
      model <- Arima(gv_ts, order = c(p, d, q), include.drift = TRUE)
  
      # Store model parameters and AIC/BIC/AICc values in matrix
      results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
  
      # Increment row index
      i <- i + 1
    }
  }
}

# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")

# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)

# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  row_spec(highlight_row, bold = TRUE, background = "#FFFF99")  # Highlight row in yellow

Comparison of ARIMA Models

p	d	q	AIC	BIC	AICc
0	1	0	2491.868	2499.012	2491.914
1	1	0	2489.464	2500.181	2489.557
2	1	0	2491.169	2505.458	2491.324
3	1	0	2493.105	2510.965	2493.338
4	1	0	2413.508	2434.941	2413.836
0	1	1	2489.862	2500.578	2489.955
1	1	1	2477.350	2491.639	2477.506
2	1	1	2475.754	2493.614	2475.987
3	1	1	2469.481	2490.914	2469.809
4	1	1	2413.555	2438.560	2413.995
0	1	2	2491.114	2505.402	2491.269
1	1	2	2461.185	2479.046	2461.419
2	1	2	2471.032	2492.465	2471.361
3	1	2	2450.991	2475.996	2451.430
4	1	2	2415.378	2443.955	2415.945
0	1	3	2479.648	2497.509	2479.881
1	1	3	2453.874	2475.306	2454.202
2	1	3	2453.964	2478.969	2454.403
3	1	3	2448.673	2477.250	2449.240
4	1	3	2416.449	2448.599	2417.161
0	1	4	2426.689	2448.122	2427.017
1	1	4	2428.569	2453.574	2429.009
2	1	4	2427.455	2456.032	2428.022
3	1	4	2429.022	2461.172	2429.734
4	1	4	2418.418	2454.140	2419.291
0	1	5	2428.602	2453.607	2429.041
1	1	5	2424.873	2453.450	2425.440
2	1	5	2426.905	2459.054	2427.616
3	1	5	2427.978	2463.700	2428.851
4	1	5	2419.521	2458.815	2420.573
0	1	6	2429.563	2458.141	2430.130
1	1	6	2426.622	2458.771	2427.333
2	1	6	2423.968	2459.690	2424.841
3	1	6	2425.868	2465.161	2426.919
4	1	6	2421.496	2464.362	2422.744
0	1	7	2430.482	2462.632	2431.194
1	1	7	2428.499	2464.220	2429.372
2	1	7	2425.761	2465.055	2426.813
3	1	7	2427.382	2470.248	2428.630
4	1	7	2423.327	2469.765	2424.789
0	1	8	2426.343	2462.065	2427.216
1	1	8	2423.969	2463.263	2425.021
2	1	8	2425.870	2468.736	2427.118
3	1	8	2420.360	2466.798	2421.822
4	1	8	2423.777	2473.787	2425.471
0	1	9	2425.562	2464.856	2426.614
1	1	9	2425.747	2468.613	2426.995
2	1	9	2428.923	2475.361	2430.385
3	1	9	2422.604	2472.615	2424.298
4	1	9	2420.926	2474.508	2422.869
0	2	0	2694.714	2698.283	2694.730
1	2	0	2589.678	2596.815	2589.724
2	2	0	2556.299	2567.004	2556.392
3	2	0	2421.927	2436.201	2422.083
4	2	0	2418.246	2436.088	2418.481
0	2	1	2483.934	2491.071	2483.980
1	2	1	2471.136	2481.841	2471.229
2	2	1	2466.091	2480.365	2466.247
3	2	1	2416.039	2433.881	2416.274
4	2	1	2409.633	2431.043	2409.962
0	2	2	2461.556	2472.261	2461.649
1	2	2	2462.294	2476.568	2462.450
2	2	2	2458.668	2476.510	2458.902
3	2	2	2416.269	2437.680	2416.599
4	2	2	2409.464	2434.442	2409.905
0	2	3	2461.613	2475.886	2461.769
1	2	3	2463.498	2481.340	2463.732
2	2	3	2448.353	2469.763	2448.683
3	2	3	2418.269	2443.248	2418.710
4	2	3	2411.237	2439.784	2411.807
0	2	4	2460.940	2478.781	2461.174
1	2	4	2444.143	2465.554	2444.473
2	2	4	2446.544	2471.523	2446.985
3	2	4	2418.292	2446.839	2418.861
4	2	4	2412.406	2444.521	2413.121
0	2	5	2421.724	2443.134	2422.053
1	2	5	2422.362	2447.340	2422.803
2	2	5	2424.094	2452.640	2424.663
3	2	5	2420.282	2452.397	2420.996
4	2	5	2414.341	2450.025	2415.218
0	2	6	2422.448	2447.427	2422.889
1	2	6	2418.950	2447.497	2419.519
2	2	6	2420.950	2453.065	2421.664
3	2	6	2420.879	2456.562	2421.755
4	2	6	2415.589	2454.841	2416.645
0	2	7	2424.378	2452.924	2424.947
1	2	7	2420.950	2453.065	2421.664
2	2	7	2423.793	2459.477	2424.670
3	2	7	2420.419	2459.671	2421.475
4	2	7	2417.589	2460.409	2418.842
0	2	8	2426.250	2458.365	2426.964
1	2	8	2422.027	2457.710	2422.903
2	2	8	2422.764	2462.015	2423.820
3	2	8	2414.724	2457.545	2415.978
4	2	8	2419.284	2465.673	2420.752
0	2	9	2422.732	2458.416	2423.609
1	2	9	2420.357	2459.609	2421.413
2	2	9	2423.369	2466.189	2424.622
3	2	9	2415.815	2462.204	2417.283
4	2	9	2419.702	2469.659	2421.402

auto.arima(gv_ts)

Series: gv_ts 
ARIMA(1,2,1)(0,0,1)[4] 

Coefficients:
          ar1      ma1    sma1
      -0.2286  -0.9327  0.4306
s.e.   0.0642   0.0314  0.0518

sigma^2 = 581.6:  log likelihood = -1205.54
AIC=2419.09   AICc=2419.24   BIC=2433.36

The parameter search suggests an ARIMA(4,2,2) while the auto.arima() suggests ARIMA(1,2,1). Lets compare these to determine which one would be better.

model_output1 <- capture.output(sarima(gv_ts,4,2,2))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1)  # Locate where coefficient details start
end_line <- length(model_output1)  # Last line of output

# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")

Coefficients: 
    Estimate     SE t.value p.value
ar1   0.0045 0.1008  0.0451  0.9641
ar2   0.0385 0.0536  0.7171  0.4740
ar3   0.0810 0.0530  1.5302  0.1272
ar4   0.5197 0.0540  9.6317  0.0000
ma1  -1.1826 0.1191 -9.9301  0.0000
ma2   0.1826 0.1182  1.5449  0.1236

sigma^2 estimated as 536.3868 on 256 degrees of freedom 
 
AIC = 9.196427  AICc = 9.197684  BIC = 9.291764 
 

model_output2 <- capture.output(sarima(gv_ts,1,2,1))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output2)  # Locate where coefficient details start
end_line <- length(model_output2)  # Last line of output

# Print the relevant section automatically
cat(model_output2[start_line:end_line], sep = "\n")

Coefficients: 
    Estimate     SE  t.value p.value
ar1  -0.2596 0.0652  -3.9825   1e-04
ma1  -0.8792 0.0360 -24.4325   0e+00

sigma^2 estimated as 708.7791 on 260 degrees of freedom 
 
AIC = 9.431818  AICc = 9.431995  BIC = 9.472677 
 

Both Residual Plots show nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.

Both Autocorrelation Functions (ACF) of the residuals reveal very few significant autocorrelations, however the plot for the ARIMA(1,2,1) does have a few more autocorrelations than the ARIMA(4,2,2).

Both Q-Q Plots indicate that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.

The Ljung-Box Test p-values for the ARIMA(4,2,2) remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals while in the ARIMA(1,2,1), all p-values remain below the 0.05 significance level implying some significant autocorrelations left in the residuals and concluding that model improvements may necessary.

Coefficient Significance: The model coefficients in the ARIMA(1,2,1) model are all significant while the coefficients for the first three AR terms and the last MA term are insignificant in the ARIMA(4,2,2) model.

fit1 <- Arima(gv_ts, order = c(4,2,2), include.drift = FALSE)
summary(fit1)

Series: gv_ts 
ARIMA(4,2,2) 

Coefficients:
         ar1     ar2    ar3     ar4      ma1     ma2
      0.0045  0.0385  0.081  0.5197  -1.1826  0.1826
s.e.  0.1008  0.0536  0.053  0.0540   0.1191  0.1182

sigma^2 = 549:  log likelihood = -1197.73
AIC=2409.46   AICc=2409.9   BIC=2434.44

Training set error measures:
                    ME     RMSE      MAE         MPE      MAPE      MASE
Training set -1.015172 23.07227 17.61682 -0.03459217 0.8088008 0.3449994
                     ACF1
Training set -0.002555783

fit2 <- Arima(gv_ts, order = c(1,2,1), include.drift = FALSE)
summary(fit2)

Series: gv_ts 
ARIMA(1,2,1) 

Coefficients:
          ar1      ma1
      -0.2596  -0.8792
s.e.   0.0652   0.0360

sigma^2 = 714.2:  log likelihood = -1232.57
AIC=2471.14   AICc=2471.23   BIC=2481.84

Training set error measures:
                    ME     RMSE      MAE        MPE      MAPE      MASE
Training set 0.2117617 26.52199 20.79503 0.01694402 0.9670959 0.4072402
                    ACF1
Training set -0.03802343

Looking at the error measurements, both the MAE and the RMSE are smaller in the ARIMA(4,2,2) model.

naive <- naive(gv_ts,h=5)
mean  <- meanf(gv_ts,h=5)
drift <- rwf(gv_ts,drift=TRUE,h=5)
accuracy(naive)

                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
Training set 8.835734 28.79736 21.69429 0.4615319 1.017407 0.4248508 -0.1290414

accuracy(mean)

                        ME     RMSE      MAE      MPE     MAPE     MASE
Training set -9.354218e-14 661.1524 594.0957 -11.0776 31.14832 11.63449
                  ACF1
Training set 0.9869592

accuracy(drift)

                        ME     RMSE      MAE        MPE      MAPE      MASE
Training set -1.085001e-13 27.40835 21.14088 0.04062096 0.9861303 0.4140131
                   ACF1
Training set -0.1290414

Furthermore, both models beat the benchmark models. Overall, it seems like the ARIMA(4,2,2) model is a better fit.

Now we can finally generate our forecast with the ARIMA(4,2,2) model.

# Generate the forecast
forecast_result1 <- forecast(fit1, h = 20)

# Plot the forecast
autoplot(forecast_result1) +
  labs(title = "ARIMA(4,2,2) Forecast",
       x = "Time",
       y = "Real Education Spending per Capita") +
  theme_minimal()

This forecast suggests that real education spending per capita will continue to increase over time.

Educational Confidence

Now lets analyze confidence in the education system, since this is our primary variable of interest.

coneduc_ts_df <- gss_num %>%
  select(year, mean_coneduc) %>%
  arrange(year) %>%
  complete(year = full_seq(year, 1)) %>% 
  mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE)) %>%
  slice(-1)

coneduc_ts <- ts(coneduc_ts_df$mean_coneduc,
                 start = min(coneduc_ts_df$year),
                 end   = max(coneduc_ts_df$year),
                 frequency = 1)

Since this data is annual, we can’t use the decompose() function for a classical decomposition. Instead, we will extract the trend and remainder parts separately.

### Trend 
trend <- stats::lowess(coneduc_ts)$y
remainder <- coneduc_ts - trend

# Plots
plot(coneduc_ts, type = "l",
     main = "Data",
     xlab = "Time", ylab = "Confidence in Education")

plot(trend, type = "l",
     main = "Trend Component",
     xlab = "Time", ylab = "Trend")

plot(remainder, type = "l",
     main = "Remainder Component",
     xlab = "Time", ylab = "Remainder")

The trend component shows a clear downward trend in sentiment and the remainder shows some fluctuations around zero with some small spikes in the early years of the survey.

gglagplot(coneduc_ts, do.lines=FALSE, lags = 9) +
  xlab("Lags") + ylab("Yt") + 
  ggtitle("Lag Plot for confidence in education") + 
  theme(axis.text.x=element_text(angle=45, hjust=1))

The lag plots show a slight relationship between consecutive observations, suggesting there may not be very much autocorrelation.

ggAcf(coneduc_ts, 50) + 
  ggtitle("Autocorrelation Function for\nAverage Confidence in the Education System") + 
  theme_minimal()

The ACF plot, however, shows some autocorrleation. Lets confirm stationarity with an augmented Dickey Fuller test.

tseries::adf.test(coneduc_ts)

    Augmented Dickey-Fuller Test

data:  coneduc_ts
Dickey-Fuller = -2.5281, Lag order = 3, p-value = 0.362
alternative hypothesis: stationary

Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our conclusions from the ACF plot.

Lets try to correct this with differencing.

diff_1 <- diff(coneduc_ts)
ggAcf(diff_1,50) +
  ggtitle("ACF of First-Order Differenced Series") +
  theme_minimal()

ggPacf(diff_1,50) +
  ggtitle("PACF of First-Order Differenced Series") +
  theme_minimal()

The ACF and PACF plots for the differenced series indicate that the data is now stationary and does not need additional differencing.

Now that we know that one difference is necessary, lets move on to a parameter search.

# Define parameter ranges
p_range <- 0:6
d_range <- 1
q_range <- 0:9  

# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)

# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)

# Initialize index for matrix row
i <- 1

# Loop through combinations of ARIMA model parameters
for (d in d_range) {
  for (q in q_range) {
    for (p in p_range) {
  
      # Fit ARIMA model with specified (p,d,q)
      model <- Arima(coneduc_ts, order = c(p, d, q), include.drift = TRUE)
  
      # Store model parameters and AIC/BIC/AICc values in matrix
      results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
  
      # Increment row index
      i <- i + 1
    }
  }
}

# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")

# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)

# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  row_spec(highlight_row, bold = TRUE, background = "#FFFF99")  # Highlight row in yellow

Comparison of ARIMA Models

p	d	q	AIC	BIC	AICc
0	1	0	-146.1814	-142.3177	-145.9314
1	1	0	-154.1756	-148.3801	-153.6650
2	1	0	-160.2939	-152.5666	-159.4243
3	1	0	-161.7735	-152.1144	-160.4402
4	1	0	-160.6396	-149.0486	-158.7305
5	1	0	-159.4146	-145.8918	-156.8099
6	1	0	-161.3818	-145.9272	-157.9532
0	1	1	-157.4889	-151.6935	-156.9783
1	1	1	-155.6315	-147.9042	-154.7619
2	1	1	-164.4317	-154.7726	-163.0984
3	1	1	-162.4820	-150.8911	-160.5729
4	1	1	-160.5157	-146.9929	-157.9110
5	1	1	-159.1041	-143.6495	-155.6755
6	1	1	-159.4157	-142.0293	-155.0255
0	1	2	-155.9474	-148.2201	-155.0778
1	1	2	-158.4396	-148.7804	-157.1062
2	1	2	-162.4920	-150.9010	-160.5829
3	1	2	-162.4265	-148.9037	-159.8218
4	1	2	-160.9624	-145.5078	-157.5338
5	1	2	-159.1660	-141.7795	-154.7757
6	1	2	-159.7447	-140.4264	-154.2447
0	1	3	-158.2586	-148.5995	-156.9253
1	1	3	-159.5933	-148.0023	-157.6842
2	1	3	-162.3082	-148.7854	-159.7035
3	1	3	-161.0628	-145.6082	-157.6342
4	1	3	-158.0646	-140.6782	-153.6744
5	1	3	-156.1666	-136.8484	-150.6666
6	1	3	-159.5780	-138.3280	-152.8088
0	1	4	-162.1637	-150.5727	-160.2546
1	1	4	-162.4235	-148.9007	-159.8188
2	1	4	-160.5475	-145.0929	-157.1189
3	1	4	-159.1997	-141.8133	-154.8095
4	1	4	-158.1940	-138.8757	-152.6940
5	1	4	-156.2203	-134.9702	-149.4510
6	1	4	-160.4167	-137.2348	-152.2062
0	1	5	-162.2605	-148.7377	-159.6559
1	1	5	-161.2509	-145.7963	-157.8224
2	1	5	-159.2674	-141.8809	-154.8771
3	1	5	-159.8130	-140.4947	-154.3130
4	1	5	-158.7977	-137.5476	-152.0284
5	1	5	-157.0933	-133.9114	-148.8828
6	1	5	-155.3629	-130.2491	-145.5250
0	1	6	-160.3046	-144.8500	-156.8760
1	1	6	-159.2570	-141.8705	-154.8667
2	1	6	-157.4700	-138.1517	-151.9700
3	1	6	-156.3756	-135.1255	-149.6063
4	1	6	-157.4528	-134.2709	-149.2422
5	1	6	-156.2534	-131.1397	-146.4156
6	1	6	-157.7621	-130.7165	-146.0954
0	1	7	-160.7412	-143.3547	-156.3509
1	1	7	-158.7809	-139.4626	-153.2809
2	1	7	-160.6976	-139.4475	-153.9284
3	1	7	-158.6980	-135.5161	-150.4875
4	1	7	-156.3615	-131.2478	-146.5237
5	1	7	-154.9472	-127.9016	-143.2805
6	1	7	-155.7651	-126.7877	-142.0508
0	1	8	-158.8138	-139.4956	-153.3138
1	1	8	-159.0816	-137.8316	-152.3124
2	1	8	-158.6981	-135.5162	-150.4875
3	1	8	-156.6980	-131.5842	-146.8601
4	1	8	-157.0020	-129.9564	-145.3353
5	1	8	-152.7827	-123.8053	-139.0684
6	1	8	-154.9958	-124.0866	-138.9958
0	1	9	-161.5029	-140.2528	-154.7337
1	1	9	-159.6613	-136.4794	-151.4508
2	1	9	-158.4304	-133.3166	-148.5925
3	1	9	-152.4592	-125.4136	-140.7925
4	1	9	-155.0022	-126.0249	-141.2880
5	1	9	-153.0854	-122.1762	-137.0854
6	1	9	-153.0750	-120.2339	-134.5295

auto.arima(coneduc_ts)

Series: coneduc_ts 
ARIMA(2,1,1) with drift 

Coefficients:
          ar1      ar2     ma1    drift
      -1.0841  -0.6679  0.6965  -0.0074
s.e.   0.1485   0.1162  0.1436   0.0038

sigma^2 = 0.002027:  log likelihood = 87.22
AIC=-164.43   AICc=-163.1   BIC=-154.77

According to the parameter search and the auto.arima() function, our best model is a simple ARIMA(2,1,1). Lets quickly look at the model diagnostics.

model_output1 <- capture.output(sarima(coneduc_ts,2,1,1))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1)  # Locate where coefficient details start
end_line <- length(model_output1)  # Last line of output

# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")

Coefficients: 
         Estimate     SE t.value p.value
ar1       -1.0841 0.1485 -7.3013  0.0000
ar2       -0.6679 0.1162 -5.7458  0.0000
ma1        0.6965 0.1436  4.8502  0.0000
constant  -0.0074 0.0038 -1.9609  0.0558

sigma^2 estimated as 0.001867907 on 47 degrees of freedom 
 
AIC = -3.224151  AICc = -3.207101  BIC = -3.034756 
 

The Residual Plot shows nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.

The Autocorrelation Function (ACF) of the residuals reveals no significant autocorrelations.

The Q-Q Plot indicates that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.

The Ljung-Box Test p-values remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals.

Coefficient Significance: The model coefficients in the ARIMA(2,1,1) model are all significant

Overall, this seems like a strong model. Lets compare it to some benchmark methods.

fit <- Arima(coneduc_ts, order = c(2,1,1), include.drift = FALSE)
summary(fit)

Series: coneduc_ts 
ARIMA(2,1,1) 

Coefficients:
          ar1      ar2     ma1
      -1.0610  -0.6314  0.7172
s.e.   0.1504   0.1203  0.1375

sigma^2 = 0.002134:  log likelihood = 85.43
AIC=-162.85   AICc=-161.98   BIC=-155.13

Training set error measures:
                      ME       RMSE        MAE        MPE     MAPE      MASE
Training set -0.01067867 0.04438487 0.03073701 -0.5315693 1.437563 0.9150038
                    ACF1
Training set -0.01935951

naive <- naive(coneduc_ts,h=5)
mean  <- meanf(coneduc_ts,h=5)
drift <- rwf(coneduc_ts,drift=TRUE,h=5)
accuracy(naive)

                       ME       RMSE        MAE        MPE     MAPE MASE
Training set -0.006753008 0.05591304 0.03359222 -0.3509282 1.552349    1
                   ACF1
Training set -0.4017722

accuracy(mean)

                        ME       RMSE        MAE        MPE     MAPE     MASE
Training set -1.238466e-16 0.09197785 0.06769079 -0.1854871 3.184123 2.015073
                  ACF1
Training set 0.7513554

accuracy(drift)

                       ME       RMSE        MAE         MPE     MAPE      MASE
Training set 3.047601e-17 0.05550374 0.03343843 -0.03171947 1.541556 0.9954219
                   ACF1
Training set -0.4017722

Our ARIMA(2,1,1) beats all the benchmark methods according to the RMSE and the MAE.

Finally, we can plot our forecast.

# Generate the forecast
forecast_result <- forecast(fit, h = 5)

# Plot the forecast
autoplot(forecast_result) +
  labs(title = "ARIMA(2,1,1) Forecast",
       x = "Time",
       y = "Average Confidence in Education") +
  theme_minimal()

The ARIMA predicts that confidence in education will remain low over time.

Conclusion

These analyses imply that that despite the government spending more on education per capita, people believe that education should be prioritized more. One possibility is that government spending on education is being misused and isn’t have as big of an effect on the quality of education as it should. Another possibility is that government spending in general hasn’t been increasing as much as it should so even if the portion of spending that is being used for education increases, the actual amount being spent isn’t meeting the demand.